NOISE-INDUCED CANARD AND MIXED-MODE OSCILLATIONS 
IN LARGE STOCHASTIC NETWORKS WITH MULTIPLE 

TIMESCALES 
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Abstract. We investigate the dynamics of large stochastic networks with different timescales 
and nonlinear mean-field interactions. After deriving the limit equations for a general class of network 
models, we apply our results to the celebrated Wilson-Cowan system with two populations with or 
without slow adaptation, paradigmatic example of nonlinear mean-field network. This system has 
the property that the dynamics of the mean of the solution exactly satisfies an ODE. This reduction 
allows to show that in the mean-field limit and in multiple populations with multiple timescales, noise 
induces canard explosions and Mixed- Mode Oscillations on the mean of the solution. This sheds new 
light on the qualitative effects of noise and sensitivity to precise noise values in large stochastic 
networks. We further investigate finite-sized networks and show that systematic differences with the 
mean-field limits arise in bistable regimes (where random switches between different attractors occur) 
or in mixed-mode oscillations, were the finite-size effects induce early jumps due to the sensitivity of 
the attractor. 
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Introduction. At functional scales, cortical networks are composed of a large 
number of cells subject to an intense noise. Such networks process and encode in- 
formation collectively, and in the variety of information transmitted, synchronous 
activity and periodic behaviors are ubiquitous [THl dHlj both in physiological (as- 
sociated to memory and attention) and pathological (epilepsy) states. Oscillatory 
patterns observed can be extremely complex. Theoretical descriptions of such com- 
plex oscillations have been satisfactorily performed in deterministic models of single 
neurons, and are generally attributed to the presence of several very different time 
scales in the dynamics. The analysis of these systems shed light on important and 
sensitive dynamical scenarios such as canard explosions (see for example |20p. The 
question we address in this work is whether such oscillatory patterns can arise in the 
context of large-scale neuronal aggregates, including a very large number of neurons, 
strongly interconnected and subject to noise. In other words, we will be interested in 
complex oscillatory patterns that can appear in large scale neuronal assemblies in the 
presence of noise. As was the case of deterministic small dimension model, we will 
specifically focus on the canard phenomenon and the related phenomenon of Mixed 
Mode Oscillations (MMOs). 

Canard phenomenon, known also as canard explosion is a transition from a small, 
Hopf type oscillation to a relaxation oscillation, occurring upon variation of a parame- 
ter. This transition has been first found in the context of the van der Pol equation [1] 
and soon after in numerous models of phenomena occurring in engineering and in 
chemical reactions [P!. A common feature of all these models is the presence of time 
scale separation (slow and fast variables). A particular feature of canard explosion is 
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that it takes place in a very small parameter interval. For the van der Pol system, 
where the ratio of the timescales is given by a small parameter e, the width of this 
parameter interval can be estimated by 0(exp(— c/e), where c > is a fixed positive 
constant. The transition occurs through a sequence of canard cycles whose striking 
feature is that they contain segments following closely a stable slow manifold and 
subsequently an unstable slow manifold. 

The work on canard explosion led to investigations of slow/fast systems in three 
dimensions, with two slow and one fast variables. In the context of these systems 
the concept of a canard solution or simply canard has been introduced, as a solution 
passing from a stable to an unstable slow manifold [31 [221 ■ Canards arise near so 
called folded singularities of which the best known is folded node [21 [211 [30] ■ Unlike in 
systems with one slow variable canards occur robustly (in 0(1) parameter regions). 
Related to canards are Mixed Mode Oscillations (MMOs), which are trajectories that 
combine small oscillations and large oscillations of relaxation type, both recurring in 
an alternating manner. Recently there has been a lot of interest in MMOs that arise 
due to the presence of canards, starting with the work of Milik, Szmolyan, Loeffelmann 
and Groeller [53]. The small oscillations arise during the passage of the trajectories 
near a folded node or a related folded singularity. The dynamics near the folded 
singularity is transient, yet recurrent: the trajectories return to the neighborhood of 
the folded singularity by way of a global return mechanism [8]. 

Due to the sensitivity of this phenomena, one may expect that canard explosions 
disappear in the presence of noise. This is implied by the results of [5], who study 
a system with small noise. Using a pathwise approach, they show that small noise 
perturbations dramatically disturb the behavior of the slow- fast system they consider, 
hiding small oscillations in the random fluctuation and having a positive probability 
to escape from the canard solution. In this picture, it seems relatively hopeless to 
obtain at the levels of noise typical of the cortical activity, and at the dimension of 
the system, collective canard phenomena and individual canards with a high degree of 
reproducibility and regularity. However, in large scale stochastic systems, averaging 
effects can take place, regularizing the probability distribution of the macroscopic 
activity, as well as that of individual cells. These macroscopic behaviors and the 
averaging effects arising at the level of large networks are mathematically characterized 
by limit systems, the mean-field limits. In the present manuscript we show that canard 
explosion and complex oscillatory pattern appear in the presence of (not necessarily 
small) noise in these mean-field limits. Moreover, we show that such phenomena may 
arise through noise induced bifurcations. 

Usual computational neurosciences descriptions of large scale networks rely on 
heuristic models which were derived in the seminal works of Wilson, Cowan and 
Amari [I] [H [3T1 [32] . These models implicitly make the assumption that noise cancels 
out through averaging effects in the asymptotic regime where the number of neurons 
goes to infinity. The analysis of these heuristic models successfully accounted for a 
number of cortical phenomena such as spatio-temporal pattern formation in spatially 
extended models [El [H [7] . 

Indeed, in order to uncover, in large scale noisy system, the presence of fine effects 
such as complex oscillatory patterns or canard explosions, one needs to very precisely 
describe the limits of the network activity as the number of neurons increases. This 
is why we will apply a recent developed approach based on the theory of interacting 
particle systems to our networks. This approach was used in a very general context in- 
cluding nonlinear models of spiking neurons, spatial extension and propagation delays 



Stochastic Networks with multiple timescales 



3 



in |28j . The resulting macroscopic states are described through extremely complex 
equations, and in general do not provide sufficient information on the dynamics of the 
global variables. However, in a simple Amari-Hopfield type of models, it was shown 
in |28| that under the condition that initial conditions are Gaussian that the equations 
can be exactly reduced to a set of deterministic differential equations governing the 
two first moments of the solutions, and these were used to analyze the dynamics of 
the macroscopic regimes [351 127] ■ This approach is precisely the one chosen in the 
present article. 

In contrast with the previous studies in this domain, we will consider here that 
the neurons interact through a nonlinear function of a mean-field term corresponding 
to the current generated at the level of each population. This model, known as the 
Wilson-Cowan model, is biologically more relevant. Moreover, in the original article 
of Wilson and Cowan [3T], the system is introduced with a slow-fast structure that, to 
the best of our knowledge, was never analyzed in the perspective of canard explosions 
and MMOs. This model has been widely studied in the context of computational 
neuroscience, and seems particularly suited for the analysis of the dynamics of large 
networks. From the computational neuroscience viewpoint, this analysis will shed 
light on the finite-dimensional Wilson-Cowan system, since incidentally the reduced 
mean-field equations will be a particular case of Wilson-Cowan system. From the 
mathematical viewpoint, the mean-field equations arising in the case of nonlinear 
global coupling were never considered in literature and lead to a number of techni- 
cal questions that will be treated in the present article. Moreover, that model has 
the advantage that the solutions of the mean-field equations have their statistical 
expectation satisfying a system of ordinary differential equations, regardless of the 
initial condition chosen, allowing for an application of bifurcation theory to analyze 
the system. The relative simplicity of the model will allow us to go very deep into the 
mathematical analysis of large scale neuronal networks in the presence of noise and 
different timescales, and address the problem of the presence of complex oscillatory 
patterns and canards in such systems. 

The paper is organized as follows. Section [T] is devoted to the introduction of the 
models and the derivation of the mean-field equations. This section deals with general 
neuron models with nonlinear synapses, multiple populations, timescales and noise. 
We apply this analysis to the Wilson-Cowan system and derive the related moment 
equations, a set of ordinary differential equations. The obtained equations are then 
analyzed in section [2j and we show that the system present canard explosions that 
can lead, in different contexts, to mixed-mode oscillations as the noise parameter is 
varied. These phenomena are recovered in numerical simulations of the finite-size 
network equations as shown in section [3] where finite-size effects are also analyzed. 

1. Network models and their mean field limits. This section introduces a 
general model of multi-population neuronal network and derives the related macro- 
scopic mean-field equations arising as a limit when the network size goes to infinity. 

Models. We are interested in the large scale behavior arising from the nonlinear 
coupling of a large number N of stochastic diffusion processes representing the mem- 
brane potential of neurons. Each neuron belongs to a specific population a € {1 ■ ■ ■ P}, 
characterizing the dynamics of neurons and their time constants, and the interaction 
with other neurons. Each population is composed of a large number of neurons de- 
noted {Na,Oi = 1- ■ ■ P). The state of neuron i in population a is described by a 
variable XI ^ E = typically describing the membrane potential of the cell and 
different ionic concentrations. The intrinsic dynamics of this variable (disregarding 
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input and interactions with other neurons) is governed by a nonhnear function fa and 
a time constant Tq. 

The action of neurons of population /3 onto neurons of population a will be mod- 
eled through a generalized version of Wilson-Cowan system, as a nonlinear sigmoidal 
transformation of a the average of the activity of all cells, weighted by the typical 
synaptic connectivity Jap- Beyond network interactions, neurons receive inputs cor- 
responding to the sum of a deterministic current /" only depending on the population 
they belong to, and stochastic current driven by independent stochastic process with 
common law (denoted with a slight abuse of notations ), that will be considered to 
be Ornstcin-Uhlenbeck processes, taking into account the synaptic filtering. In their 
most general form, the class of models analyzed here satisfy the network equations: 




Uxi) + Sa XI V ^ V xi + /r + Ct\- (1.1) 



where Sa{x, y) represents the effect of an input current y onto a cell of population a 
with potential x, X]^=i 7^ Sp(j")=0 X} the total current received through the net- 
work, and p : N I— !■ { 1 • • • P} is the population function associating to a neuron j the 
population /3 it belongs to. 

This class of models is particularly interesting since it allows dealing with a num- 
ber of situations arising in neurosciences. A prominent example of this class that will 
be useful in our developments is the Wilson-Cowan model, or activity-based firing 
rate model [3T1 [32] ■ This model describes the mean firing-rate variable of a neuron (a 
real variable) of neuron i in the population a through the differential equation: 



p 



0=1 p 



Xl+I?+S,l) (1.2) 



/3=1 ^ p(i)=/3 

In this model, the variable XI represents the activity of the cell, and the right hand 
side is the incoming firing-rate to the cell i, filtered through the voltage to rate function 
Sa characterizing the integration properties of the synapse when subject to a input 
current (see [l3]). 

Another model of this class allows taking into account finer phenomena arising 
in the cell's activity motivating the description of additional variables, such as slow 
adaptation currents [11' or slow modulation currents |17j . These models correspond to 
minimal descriptions of systems with slow adaptation proposed by Izhikevich in |17j : 
in two-dimensional systems, a variety of behavior can be accounted for using a single 
slow adaptation. We slightly generalize this model and define a multi-population 
networks with mean- field slow modulation current. These systems are described by 
the following equations, for neuron i in population a: 



St' - + Sa{Y.0 = l Y.p(j)=0 Xl + XalJl + If + Ct^) 

I =ea{k + ^ui - 1^ XI) 



Beyond these simple models that will be exactly reducible to ODEs, the class 
of models includes classical models in neurosciences such as the Fitzhugh-Nagumo 
model [X e M^, has a cubic and an affine component, and S{x, y) is either a linear 
function when considering electrical synapses or a sigmoid transform when considering 
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chemical synapses) or the Hodgkin-Huxley system. These models differ from usual 
systems analyzed in the mean-field limit in that there is a nonlinear summation of 
the interactions, here through our voltage-to-rate function. Limits of such interacting 
stochastic processes were the number of neurons tends to infinity were never addressed 
for such interactions considering the nonlinear transform S. We now derive the limits 
of these stochastic networks equations as the number of neurons go to infinity. 



Macroscopic Limits. Based on the general network given by equations (1.1), we 
derive the macroscopic mean-field limit under a few assumptions on the dynamics. For 
the sake of simphcity, we denote by /(x) : (resp. S{x, y) : E^ x E^ H' E^) 

the function with component a equal to fa{x") (resp. Sa{x°' ,y°'))- We make the 
following assumptions: 
(HI). / is X-Lipschitz continuous 

(H2). / satisfies the linear growth condition: f{x) < C(l -|- |a;p). 
(H3). S is L-Lipschitz continuous in both variables 
(H4). S is bounded with supremum denoted ||S'||oo- 
Remarks. 

• Rigorously, in the Fitzhugh-Nagumo or Hodgkin-Huxley system, f is not globally 
Lipschitz continuous and does not satisfies a linear growth condition. However, these 
functions are locally Lipschitz continuous, and satisfy the estimate x f{x) < 7^(1-1- 
|a;|^). These two assumptions allow to demonstrate the same kind of properties, 
by using ltd formula and truncation method. This was done in J28^ in a more 
general context, and the interested reader may readily apply this generalization to 
our system. 

• In these assumptions, we consider without loss of generality the same constants K 
for the sake of simplicity. 

• All these assumptions are obviously satisfied for the Wilson-Cowan system. 

We will show that in the mean-field limit, the propagation of chaos applies and 
that the state of all neurons in that limit satisfy the implicit equation: 

d X" I ^ r n 

which can be written in vector form: 

^ - fiXt) + S {Xt, J ■ E [Xt] +It+ 6) (1-5) 

where J = {Jap)c.,pe{i-P} e M^^^ and It = (/r)a=i...p, 6 = (Cf , a e {1 ■ ■ ■ P}). 

This equation is not a usual stochastic differential equation in that in involves 
the expectation of the solution. It is rather an implicit equation on the space of 
probability distributions. The first question that may arise in this type of unusual 
stochastic equation is its the well-posedness. This is subject to the following: 

Proposition 1.1. There exists a unique, square integrable, strong solution to 



the equations (1.4) starting from a square integrable initial condition. 

Proof. The existence and uniqueness of solutions is demonstrated using a routine 
contraction argument. Let us fix a square integrable initial condition random variable 
Xq G E^ and define the map $ acting on the space of square integrable stochastic 
processes M^{[0,T], E^): 



: X ^ (^Xo + (/(X,) + S{Xs, J ■ E [Xs] +Is+ is)) ds 



te[o,T] 
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Fixed points of $ are exactly the solutions of the mean-field equations. Moreover, for 
X a square integrable, we have ^{X) G M'^i[0,T], E^). Indeed, thanks to Cauchy- 
Schwarz inequality, we have: 



E 



sup mx)t 

te[o,T] 



< 



4|e[|XoP] + T £ C'{l + E[\Xs\']) + \\S\\lds^ 



Fixing a square integrable process e A^^([0, T], E^), we define the sequence of 
square integrable processes = We note D}' = E supt^^a,T] l^t'^^ " ^fP] 

II J|| the operator norm of J and = J • E [X^] + Ig + Cs- We have: 



< ST e| sup /' \f{x:) - fix:-')\' + \s{x:, z-) - s{x:-\ z:)\' 

[ se[o,t] Jo 



+ \s{x:-\z:)-s{x:-\z:-')\uu 



<3T]S\ sup riK^ + L')\x:~x:-'\^+L%j-E[x:-x:-']\^du\ 

[ se[o,t] Jo J 
<3rE| sup r{K^ + L')\X:-X:-'\^+L^J\\E[\X:-X:-'\^] du\ 

[ sG[0,t] Jo ) 

<3T(K'^ + L\l + \\J\\)) f Dl'-^s 

Jo 

We hence obtain by immediate recursion: 

„ < (^^0 
n! 

within' = 3T{K^+L'^{1+\\J\\)). Dj. is a finite quantity because of the boundedness of 
E From this inequality, routine methods using the Bienayme-Chebyshev 

inequality allow concluding on the existence and uniqueness of solutions of the mean- 
field equation (see e.g. HSj). □ 

The next theorem proves that the unique solution to the mean-field equations will 
describe exactly the behavior of the network when the network size goes to infinity 
(i.e. when all Na — > oo), and that all neurons will be independent provided that 
their initial condition is. This latter property is referred as the propagation of chaos 
property. 



Theorem 1.2. If the initial conditions to the network equations (1.1) are inde- 



pendent, and identically distributed population per population, then in the mean-field 
limit, the propagation of chaos applies and the laws of the network solutions converge 



towards that of the mean-field equations (1.4), in the sense that for any i G N there 



exists XI with law given by the solution of the mean-field equations and such that: 



E 



sup \Xl-Xl\ 
te[o,T] 



< 



cm 



Proof. The proof of the theorem uses the usual coupling argument popularized by 
Sznitman 26] consisting in identifying an almost sure limit for the network equations. 
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In details, considering neuron i of population a, we define the coupled process XI 
solution of equation (1.4) with initial condition identical to the initial condition 
of neuron i in the finite network and driven by the same noise Q: 



dt 



^Aq — Aq 



Uxi) + [xi E^=i Jc^p^ [^f ] + It + 



We denote by J the maximal value of the coupling strengths J^^. We aim at show- 
ing that X^ almost surely converges towards X*, and to this purpose analyze the 



distance between the two processes, Dl 



E 



Because of in- 



terchangeability in the network, that quantity does not depend on the specific neuron 
i we consider but only on the neuronal population a it belongs to. For the sake of 
simplicity, we note: 



yN,a _ V^P Jag y] 



7N,I 



13=1 Si SpO )=/3 E 



In order to analyze this distance, we decompose it into the sum of four simpler terms 
that are more easily controllable: 



\xi-xi\< 



\Uxi) - uxi)\ + \s{xi, zr') - s{xi, zr')\ 

N\\ I \cilvi ^a,N 



\Xl - XI 



ljy: 



/3=1 



1 

iv7 



13 



+ \s{x:,zf^'' 

<{K + L) 

readily yielding: 

t p p 

Dt< [ {K + L)D^ + LJY^dI + LJY,^ 

a -i ail 



S{Xl,zr')\ + \S{Xl, Z,"'^) - S{Xl, z':)\ ds 
p 




E 



13=1 



13=1 



■E 



Xi-E[Xi] 

pU)=P 



ds 



and hence for Dt ~ X)q=i ^? obtain: 



Dt< {K + L + PLJ)Ds + PLj ^—1 

Jo a-^ P 



J2 xi-K[xi\ 

pU)=0 



ds (1.6) 



Let us now evaluate the second term of the righthand side. We have, using Cauchy- 
Schwarz inequality: 



E 



Y,Xl-E[Xi] 

jei3 



< E 



J2xi-^[xi\ 

je/3 



1/2 



The processes {X^ — E [X^^ ) are centered and independent since they are driven by 
independent processes and ^j. Developing the square of the sum, we are hence 
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left with less that Nj3 non-zero terms which are equal the variance of , na mely 
E \{Xi — E which is a bounded quantity as a result of proposition ll.ll Let 

us denote by C this bound. Returning to equation (1.6 1, we obtain: 



Dt< f iK + L + PLJ)D,ds + TP'^LjJ 

Jo V 



and Gronwall's lemma implies that: 

Dt < 



min^ y/Np 



oK2t 



with K2 = {K + L + P LJ) and Ki ~ P LJ\/CT . This quantity tends to zero when 
all Na tend to zero (i.e. when min^ — > 00), which implies the almost sure conver- 
gence {XI) towards {XI), and hence convergence in law of {XI) towards the solution 
of the mean-field equations. Moreover, considering a finite number I of neurons in the 
network (ii, • • • it is easy to see using the above result that {X^^ , ■ ■ ■ ,Xl') con- 
verge almost surely towards {Xl^ , XI' ) which is a set of independent processes. 
We therefore have propagation of chaos, i.e. that in the limit where all Np tend to 
infinity, the state of these neurons are independent for all times. □ 
Remark. Note that this proof has been performed in the absence of Brownian additive 
noise in contrast with usual approaches. Yet, the noisy input received by each cell allows to 
recover the same kind of properties. We also emphasize the fact that the proof provided here 
readily extends to networks with additive noise. 

As announced, the results of proposition 1 1 . 1| and theorem |1 . 2| readily apply to the 
Wilson-Cowan model, as we now show: 

Corollary 1.3. In the mean-field limit where all population sizes tend to 



infinity, the Wilson-Cowan system (1.2) in the mean-field limit satisfies the propaga- 



tion of chaos property and the law of neurons of population a asymptotically satisfy 
the mean-field equations: 



dxi 

' dt 



where we denoted fJ.p{t) = E 
equation: 



-x^ 



s 



p 

E 



Jo.p^lp{t) + If 



(1.7) 



X 



These quantities (/i^(t))^=i...p satisfy a closed 
p 

(t) + G„(^ J„^^^(t)-Mr) (1-8) 



/3=1 



where 



Go. :^x^E[So.{x + C?)]- 



In the case of slow modulation currents (1.3 I, the mean-field equations read: 



= + 5(eLi Jo.0iic.{t) + + /r + ) 



dUt 
dt 



= £a{k + -fUt - Y.B=1 M/3(i)) 



(1.9) 



(1.10) 
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and in that case, Ut is a deterministic variable satisfying an ODE and again, the 
expectation /ict(i) of the process X" satisfies a closed set of ordinary differential equa- 
tions: 



dUt 



Proof. Both cases are can be readily set in the general framework of the current 
section, hence proposition |1.1| and theorem |1.2| apply, ensuring propagation of chaos 
and yielding the announced mean field limit. The averaged equation is trivial to derive 
taking the expectation of the process. □ 

We therefore proved that the systems under consideration here enjoy a consider- 
able dimension reduction: the original system is a very large network of nonlinearly 
interacting stochastic processes with random currents, in which the mean is not so- 
lution of an ordinary differential equation, in contrast to the mean-field equations 
describing the infinite-size limit which gives rise to a closed-form ODE for the mean. 

The global dynamics of the system can hence be very precisely understood by 
analyzing the dynamics of the obtained system of differential equations. These involve 
an unspecified function, Ga, which is the expectation of our sigmoidal transform 
Sa{x + with respect to In detail, if has a density pf , this function writes: 

G^{x)= [ S{x + y)dp^{y), 

coupling the dynamics of the mean equation to the distribution of the noisy current 
received by each cell. In what follows, we will consider that the noisy current is a 
stationary Ornstein-Uhlenbeck process, i.e. are centered, Gaussian processes with 
constant standard deviation denoted CTq. In the case where Sa{x) — eri^gax), a 
simple change of variables and integration by parts (see e.g. [53]) gives the function 
Ga in closed form: 

G„(x) = erf(^=i^^). (1.12) 

The noise level is hence a parameter of the vector field, and may lead to noise- 
driven bifurcations. 

From these fully determined systems, we are now in a position to use classical 
tools from the analysis of slow-fast systems to uncover the presence of canards and 
mixed- mode oscillations in macroscopic activity in relationship to the noise level cTq. 

2. Canards and MMO in the mean-field equations. We now analyze the 
limit systems (1.8) and ( |1.11[ ) with a particular focus on the effects of noise on the 



dynamics. We start by showing the role of noise in the emergence of canard explosion, 
prior to extend the model to slow modulation currents systems where the additional 
variable gives rise to noise-induced MMOs. 

2.1. A two-populations Wilson-Cowan system: noise-induced canard 



explosion. In this section, we study the mean-field equations (1.8) arising from 
Wilson-Cowan system with two populations. The network equations are given by ( |1.2[ ),| 
with P = 2 populations evolving on vastly different timescales Ti = e and T2 = 1, 
= ze, ~ 0. Moreover, in order to reduce the system to closed-form mean equa- 
tions, we choose Sa to be erf functions and ^" to be a stationary Ornstein-Uhlenbeck 
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process with standard deviation ctq . As a result of the previous section, we know that 



the associated mean- field equations (1.81 have solutions whose mean satisfy a closed 
set of ODEs given by the equations: 

Gi{JiiUi + J12U2 + Ze) ^22) 
G2{J2lUl + J22U2) 

In these equations, the two different effective nonlinearities Gi and G2 are of the 




form (1.12) and therefore depend on the noise intensities (Tq.. 

Let us now analyze the dynamics of these equations. We start by analyzing the 
behavior of the solutions upon variation of the parameter z^^ corresponding to fixed 
applied current to the first population. The behavior of the system as a function of 
current levels Ze is particularly important in our context for its relationship with the 
slow modulation currents as will be made explicit in the following section. Moreover, 
this parameter plays the important role of moving the system from the excitable 
to the spiking regime. It is well known that such a parameter transition, under 
some degeneracy assumptions, can lead to a canard explosion, i.e. a transition from 
small amplitude oscillations to relaxation oscillations, see for example [H]. And in 
the present system, canard explosions can be found by varying the parameter z^ 



are observed (see bifurcation diagram generated with AUTO in Fig. 2.1(a)). The 
diagrams suggest that the shown canard explosions are non-degenerate (see (zlj for a 
list of non-degeneracy conditions). It would be straightforward but tedious to verify 
these conditions and we omit such computations in this paper. This explosion informs 
us of the presence of canard explosions, upon variation of the input, in the Wilson- 



Cowan system (2.1), in the presence of noise. 

The question that arises is how this explosion depends on the amplitude of the 
noise neurons are submitted to. This can be addressed here using usual bifurcation 
and singular perturbation methods since the noise appears as a parameter of the 
ODEs, acting on the maximal slope of the sigmoids. Numerical simulation of the 
bifurcation diagram shows indeed that the presence of the canard explosion depends 
on the level of noise, and in particular that noise can induce canard explosions. As 



an example of such noise induced canard explosion^ we provide in Fig. 2.1(b) the 



AUTO bifurcation diagram of equations (2.1) upon variation of cti the noise received 



by neurons in the fast population. Therefore, noise can also bring the system from 
the region of excitability to the region of spiking, through a non-degenerate canard 
explosion. 

As a side result, this analysis confirms the direct relationship between levels of 
noise and synchronized oscillations in the system [291 127) . This is an extremely sur- 
prising phenomena for which no microscopic interpretation has been provided thus 
far, and which in the context of slow-fast Wilson-Cowan system under consideration, 
appears even more surprising: indeed, the presence of this explosion shows that very 
small changes in the value of the variance of the noise yield dramatic qualitative 
switches on the macroscopic dynamics. Moreover, since we have proved that all neu- 
rons have the same probability distribution population-wise, the finding shows that 
the mean presents a sudden switch to a large-amplitude periodic orbit, correspond- 
ing to all neurons presenting an extreme change upon very small variation of the 
noise. More than noise-induced oscillations, the phenomenon shows how sensitive the 
dependence upon noise levels is in large-scale systems. 

Now that the two-population system has been analyzed as a function of parame- 
ters Ze and (7i, we proceed by adding to the same system a slow adaptation current. 
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Figure 2.1. Canard explosion in the 2D WC model, for a variation of parameter Ze (top 
panels) and v\ (bottom panels). Shown are the bifurcation branch with the two explosions (left 
panels), cycles corresponding to the canard explosion highlighted by a dashed box on the bifurcation 
diagram (right panels). 



2.2. The slow modulation current Wilson-Cowan system. We now en- 
rich the network model analyzed in the previous section by considering slow adapta- 
tion currents. The underlying microscopic equation correspond to a two-populations 



Wilson-Cowan system with slow adaptation as introduced in equation (1.3). In order 
to simplify the model, we fix Ai = 1, A2 = and = P = 0. The results of sectionjl] 
ensure that the propagation of chaos occur and that the network converges towards 



mean-field equations, given by equation (1.10), and that the mean of the solution to 
the mean-field equation satisfies an equation of type ( |1.11[ ), which in our two pop- 
ulations setting with one fast and one slow population and parameters provided in 
section 12. 11 read: 



U2 = -U2 + G2{J2lUl + J22U2) 



(2.2) 



where the slow adaptation current (noted U in the general equation ( |1.11[ )) is denoted 
by Zg. This notation clearly shows that equation (2.2) can also be seen as general- 



ization of the two population Wilson-Cowan system (2.1) obtained by turning the 



parameter into a dynamic variable. In these equations, the new parameter k con- 



trols the position of the equilibria of (2.2). One of the most common routes to MMOs 



is the passage through a so-called Folded Saddle Node of type II (FSN II, see [22]). 
This transition is similar to canard explosion in the sense that it corresponds to a 
passage from a regime of a stable excitable equilibrium to MMOs or relaxation oscil- 



lators. We now investigate the presence of FSN II in our system (2.2) by considering 
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Figure 2.2. The top left panel shows a two-parameter bifurcation diagram of the 3D Wilson 
Cowan system ( |2.2[ l in the (k,<Ti)-plane, with Hopf (black) and period-doubling (blue) bifurcation 
branches. The red curve is the locus of FSN II bifurcation. All other panels show the long-term 
behavior of the system for particular values of k with fixed ai = (panels (a) to (f)) and for 
particular values of a\ with fixed k = —1 (panels A to F), respectively. 



the reduced system, obtained by setting e = 0: 



— —Ui + Gi(JiiUi + J12U2 + Ze) 
U2 = -tU2 + G2(J2lWl + J22U2) 

z'e = k + — lii — U2 



(2.3) 
(2.4) 
(2.5) 



The constraint (2.3 1 is ehminated by solving this algebraic equation in U2- This 
provides an expression of U2 as a function of [ui, Ze): 



U2 



J12 



and Gi ^ is expressed in terms of the inverse of the erf function and of ci , dependence 
that is made explicit in F. We can now rewrite (2.3)-(2.5) in the form: 



dF 
dui 



Ui 



-F{ui,Ze,<Ti 

dF 



G2( J21U1 + J22F{ui, Ze, CTi)) 



{k + JZe -Ui - F{ui,Ze,ai)) 
Ui - F{ui,Ze,(Ti) 



Now the fold line for system (2.6l-(2.7) corresponds to the equation 



dF 
dui 



(2.6) 
(2.7) 

~ and 



folded singularities correspond to points on the fold line for which the RHS of (2.61 
is 0. An FSN II point is a point on the fold line for which the RHS of both (2.61 
and (2.7) are zero. Such points are of codimension 1 in the parameter space, and, 
under some non-degeneracy conditions, correspond to a passage of a true equilibrium 
of (2.6l-(2.7) through the fold line As a parameter goes through an FSN II point 
the system passes from a regime of excitability (stable equilibrium near the fold) to a 
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Figure 2.3. Panel (a) shows the critical manifold of system ^2.2\ (green surface) viewed in 
the three-dimensional phase space (ui,y,Ze). Also shown are both the fold curves (black lines) 
and an MMO of type 1^ , solution to the problem, whose existence can be explained by the presence 
of a folded node pfn (black dot) in the system. Panel (b) shows the time profile of ui for this MMO. 



regime of MMOs or relaxatiorj^ depending on the global part of the flow [8 . On the 
level of folded singularities, FSN II corresponds to a passage from folded saddle to 
folded node, and folded nodes are associated to the occurrence of MMOs (HUl HI El]- 
Note eventually, as demonstrated in [16 , the presence of a curve of Hopf bifurcations 
in the vicinity of an FSN II point. 

This is precisely the scenario of transition occurring in out system. Indeed, in 



Figure 2.2 we display in the upper left corner the FSN II curve (shown in red), 
computed using MATLAB based on the above formulae. Nearby (in black) is the curve 
of Hopf bifurcations (HB) for e = 0.005, which, as expected, is a close approximation 
of the FSN II curve. 

Fine analysis of the passage from excitability to MMOs shows a typical rich 
bifurcation structure |121 116j . This structure, occurring as a result of changing the 



parameter k, is highlighted in figure 2.4 In particular, we find a cascade of period- 
doubling bifurcations (PD) (best seen on panel (b), which is an enlargement of panel 
(a)) and a family of complicated closed branches or isolas (we computed four of them). 
This structure is compatible with the analysis done in [16 on a minimal system of a 
singular Hopf bifurcation. 

Changes in the parameters {k, cti) correspond to paths in this diagrams governing 



the bifurcation scenarios. Two particular paths are displayed in Figure 2.2 (i) varying 



k while keeping cti = 2 fixed (small letters) and (ii) fixed k = —1 and CTi varying 



(capital letters), and a phase space view of an MMO trajectory is given in Figure 2.3 



Multiple MMO traces corresponding to these two variations of parameters highly the 
typical traces found inside the closed curve corresponding to FSN II. The scenario (i) 
will be explored numerically in section |3] in the network equations, and scenario (ii) 
highlights the fact that for k fixed and a increasing from to a positive value, we 
observe a transition from excitability to MMOs as a noise induced bifurcation. 

Therefore, the analysis of the system of ODEs arising from taking the average of 
the solution to the mean-field equations in Wilson-Cowan systems with slow adapta- 
tion currents yields to the conclusion that (i) qualitative behaviors sensitively depend 
upon noise levels, as we already learnt from the analysis of section 2.1 and (ii) that 
very delicate phenomena such as canard-induced MMOs persist in the presence of 



^In fact other types of dynamics could occur, but a passage to MMOs or relaxation is what is 
often observed. 
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Figure 2.4. Bifurcation diagram of the original system \2.2\ as a function of parameter k. 
Panel (a) shows a large view of the branch of equilibria (black) and the branches of periodic solutions 
(blue) born at two Hopf bifurcations (black dot) that are symmetric with respect to k = 0. Panel (b) 
is an enlarged view in the parameter region where additional branches of periodic solutions exist, 
both branches born at multiple period- doubling bifurcations (PD) as well as closed branches (isolas). 
The vertical axis correspond to the computed L2-norm of the solutions. 



possibly large noise. This phenomenon is, again, very surprising, and may only be 
understood in the context of collective effects in very large networks. Indeed, large 
noise of course will destroy the fine structure of small oscillations in regimes where 
the noiseless system presents mixed-mode oscillations, but we observe here that these 
persist for large noise as a collective effect. Even more interesting is that noise ac- 
tually induces the apparition of such solutions, and that the shape of the solution 
sensitively depends on the precise noise levels neurons are submitted to. 

3. Back to the network dynamics. Thus far, we showed that in the limit 



N ^ oo, the network equations (1.2) can be reduced to mean-field equations whose 
mean satisfy an ordinary differential equation where the noise parameter, appearing 
as a parameter of the limit equation, is related to the presence of canards and MMOs. 



We now return to the original network equations (1.2 1 with finite N. We show that 
the predicted presence of complex oscillations as well as canard explosion are generally 
recovered in the simulations of large but finite networks. In certain parameter regions, 
systematic errors related to the finiteness of the number of neurons considered however 
appear. These finite-size effects, corresponding to the presence of random MMOs in 
two-populations networks (in which cases no MMO is present in the mean-field limit) 
are analyzed in view of the competition between multiple attractors of the mean-field 
limit and finite-size fluctuations around that limit. In three populations networks, the 
MMO solutions are generally very fragile close from the transitions, and the finiteness 
of the number of neurons generally induces early jumps shortening the duration of 
the small oscillations phase. 

The simulation of large networks presents a number of theoretical and technical 
obstructions. Indeed, at least three difficulties arise: the stiffness of the solutions 
and the large periods of relaxation oscillations related to the slow-fast nature of the 
equations, the dimensionality issue related to the number of neurons considered (iV), 
and the stochastic nature of the trajectories. The stiffness issue is generally taken 
care by using refined simulation algorithms, that generally do not have a stochastic 
counterpart. Here, accurate simulations will necessitate to use small timesteps dt 
and large simulation times T (due to the typical periods of relaxation oscillations). 
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which will raise two difficulties: the computation time will increase proportionally to 
T/dt, and the number of independent random variables we need to draw during the 
simulation will increase proportionally to NT/dt, which can raise concerns about the 
fact that the correlations of the pseudo-random sequence of random numbers may 
surface. Moreover, with large simulation times, error accumulation will arise, pro- 
portionally to the computation time. With these limitations in mind, we developed 
a vectorized simulation algorithm based on Euler-Maruyama scheme implemented in 
Matlab@, which remains relatively efficient (computation times are almost indepen- 
dent of the network size) at the levels of precision needed. Typically, the simulation of 
a N = 2 000 neurons network over a time interval T = 1 500 with timestep dt = 0.01 
will take 66s on a Macbook Pro with processor 2.3 GHz intel core 17 with 16GB 1600 
MHz DDR3. The specific development of computational algorithm, precise analysis 
of the error ad limitations raised by the number of random numbers drawn are not 
in the scope of the present paper, and we are confident that our implementation, 
with proper choices of timestep, total times and number of neurons, provide accurate 
descriptions of the network behavior in the examples provided below. 

3.1. Canard Explosion in the 2 populations Wilson-Cowan system. We 

start by analyzing the finite-size effects in our two population network. For our choice 
of parameters, we have seen that a subcritical Hopf bifurcation and canard explosion 



occurs when varying the intensity of the noise (see Fig. 3.1(a) I. The full bifurcation 



diagram, plotted in Fig. 2.1 displayed a canard explosion with a branch of limit cycles 
persisting until very large values of the noise standard deviations (above 45). Such 
standard deviations correspond to extremely large noise that are generally very hard 
to simulate accurately and seem unreasonably large in the Wilson-Cowan model. A 
restricted bifurcation diagram is displayed in Fig. 3.1(a)| and shows that, as a function 



of noise levels, the system either presents a unique stable stationary distribution (blue 
region), an unstable fixed point and a stable periodic orbit corresponding to relaxation 
oscillations for large noise (pink region), and an intermediate noise region (yellow) 
where we have co-existence of a stable stationary distribution and a stable periodic 
orbit, separated by an unstable periodic orbit. This parameter region is referred to 
as the bistability regime. 

In the case where a unique stable solutions exists for the limit averaged equa- 
tion (pink or blue regions), network trajectories will display random fluctuations well 
approximating the attractor. This is what is plotted in Fig. |3.1[ both in the case 
of stationary and periodic solutions. The canard explosion, displayed in section [2] 
corresponds to the sharp transition from fixed point to relaxation oscillations, and 
due to the fact that the canard explosion occurs through subcritical Hopf bifurcation 
prevents from observing the full canard explosion and canard cycles, that correspond 
to unstable trajectories of the limit ODE and of the network equations. 

In the bistability regime, trajectories are more complex and finite-size effects will 
result in random switches between the two attractors. Typical shape of the phase 



plane, displayed in Fig. 3.2(a) show the coexistence of relaxation oscillations and 
stable fixed points. The separatrix between the two attractors is the unstable limit 
cycle arising from the Hopf bifurcation, winding around the fixed point. Incidentally, 
we observe that the unstable periodic orbit and the relaxation oscillations cycles are 
extremely close (they are arbitrarily close when the intensity of the noise approaches 
the value corresponding to the fold of limit cycles). Therefore, in this parameter re- 
gion, random switches from relaxation oscillations to random fluctuations around the 
stable fixed point will occur, phenomenon that we term random attractor switching. 
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(a) Bifurcation Diagram (b) Stationary regime, v = 0.5 




(c) Periodic Regime, v = 2 



Figure 3.1. Network Simu latio ns for the Wilson-Cowan two populations system, whose asymp- 
totic dynamics was analyzed in \2.1\ (a) zoomed bifurcation diagram for relevant noise intensities, 
(b) and (c) are examples of trajectories of the network equations in the cyan and pink regions re- 
spectively, (b): stationary solutions, v = 0.5; (left) 10 randomly chosen trajectories of the network 
(blue: population 1, red: population 2) display random fluctuations of the around the predicted 
stationary average (cyan: population 1 and black: population 2) computed through the limit ODE. 
(right) empirical average vs theoretical limit average, show a very precise fit. (c) Same simula- 
tion and color code as in (b) but for v = 2 corresponding to the oscillatory regime (pink region). 
Simulation parameters: 2 000 neurons per population, Euler-Maruyama scheme with dt = 0.01 and 
randomly chosen initial conditions. 



Similarly to what was shown in [6], we may interpret the switching trajectories 
as random MMOs. However, in contrast to their setting, the noise perturbation 
involved in our case are not driven by white noise with constant standard deviation. 
As we analyze in a forthcoming study, the finite-size perturbations correspond to 
Gaussian perturbations with standard deviation function of the solution of the mean- 
field equations. The nature of this process implies that the random switches do 
not occur completely randomly, but depend on the size of the fluctuations, which is 
governed by a function of the mean-field equation. A consequence will be, as presented 
in the example plotted in Fig. 3.2(b)[ that the system has a bias to switch in phase 
with the virtual periodic attractor. This surprising phenomenon will be analyzed in 
greater detail in a forthcoming study. 

Let us eventually discuss the occurrence of switches as a function of the size of the 
network and of the parameters. We observe that the attraction bassin of the stationary 
solution becomes increasingly large as noise approaches the value corresponding to the 
fold of limit cycles, and increasingly small when approaching the Hopf bifurcation. 
Close from the value of the fold of limit cycles, random switches from the periodic to 
the stationary solutions will become very probable and switches from the stationary to 
the periodic solutions are less probable, implying that the averaged time spent in the 
neighborhood of the stationary solution will become larger as we approach the fold of 
limit cycle. A symmetrical phenomenon will arise close from the Hopf bifurcation. In 
order to illustrate this phenomenon, we display in Fig. 3.3 the average time spent in a 
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0.5 1 1.5 2 

(a) Phase plane, v = 1.1 




(b) Trajectories, v = 1.1 



Figure 3.2. Bistable regime in the two-populations Wilson-Cowan system, (a) Geometry of the 
phase plane of the limit ODE, shows the presence of a relaxation oscillations cycle (black) and an 
unstable cycle (blue) delineating the attraction basin of the stable fixed point. Red and Green curves 
are the nullclines of the system, the triangle (square) depicts a saddle (unstable focus) fixed point, 
(b) depicts a typical trajectory of random attractor switching, with 1000 neurons per population and 
V = 1.1, same color code as in Fig. \3. 1\ Dotted lines depict the location of the stationary mean for 
both populations. 



neighborhood of the stationary solution (actually inside the deterministic attraction 
bassin of the stable fixed point depicted in blue in Fig. 3.2(a) ) as a function of the 
noise level, averaged across 20 trajectories. These switching probabilities also depend 
on the size of the network. Indeed, as the number of neurons tends to infinity, that 
probability tends to zero, and the trajectory will be deterministically locked to that 
predicted by the mean-field theory and the initial condition, and the smaller the 
network size, the more likely the switches will occur. 



Figure 3.3. Proportion of time along trajectories spent in the vicinity of the stationary solution 
as a function of noise levels. Colors correspond to the regions in the bifurcation diagram \3.1(a)\ 



We therefore conclude from the present analysis that finite-size effects barely 
affect the solution, and that the mean-field model accurately describes the behavior 
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of the network in most of the parameter space, except during the canard explosion, 
in an exponentially small region of parameters, and in that region, random switches 
between multiple attractors may occur, as a function of the network size and the 
geometry of the limit dynamical system. 

3.2. Mixed-Mode Oscillations: early jumps and hidden small oscilla- 
tions. We now investigate the behavior of finite-sized networks in the 2 populations 
Wilson-Cowan system with slow adaptation currents analyzed in the mean-field limit 
in section [2?2| The analysis the resulting mean-field system showed that the system 
presented a FSN II bifurcation related to the presence of mixed-mode oscillations, the 
shape of which sensitively depended of the noise levels cti and k, as depicted in Fig. |2.2[ 
The reduction obtained of these equations allowed analyze the behavior for extremely 
high values of the noise which are not accessible through numerical simulations of the 
noisy network, chiefly due to precision constraints as we mentioned. 

We now confront the different behaviors observed in section 12.21 to simulations 
of the noisy network equations with the same parameters. We expect in this case 
that the trajectories obtained in the MMO regime to be less resilient to perturbations 
than the stationary and periodic orbits analyzed in the two-populations network, 
due to the sensitivity of the attractor in particular in the small oscillations phase. 
In contrast with the analysis of the two-populations Wilson-Cowan model, we will 
therefore observe more drastic finite-size effects even in regimes where there exists a 
unique attractor, since it is the very structure of the attractor that presents sensitivity 
to perturbations. These will lead to escape the attractor during the small oscillations 
phase, or early jumps. 

This is what we observe in Fig. |3.4| where we reproduce the bifurcation scenario 
displayed in Fig. |2.2| for fixed noise ai = 2 and different values of the parameter k 
(the only scenario we did not display was (B) , which requires very simulations on very 
long time intervals due to the very large period of the mixed-mode oscillation) . In all 
cases, we recover the qualitative nature of the solutions predicted in the infinite-size 
limit, and early jumps, corresponding to systematic random escape of the attractor, 
occur due to finite-size effects. Moreover, we observe that the small oscillations tend 
to be hidden by the noise related to the finite-size effects. These two finite-size 
effects of course disappear when considering larger networks and one recovers with 



precision the attractor predicted by the theory, as we show in Fig. 3.5 for a network of 
N = 10 000 neurons. But as soon as the number of neurons is smaller, the network will 
systematically fail to follow precisely the predicted MMO and will generally escape 
the attractor prior than the infinite-sized system. 

4. Conclusion. In a number of applied domains such as physics, chemistry and 
neurosciences, the emergence of complex oscillations such as mixed-mode oscillations 
have been theoretically related to canard explosions and multiple timescales dynamics 
of deterministic dynamical systems (see 1121 for a review). But the mechanisms un- 
derlying the existence of such orbits generally involve a number of interacting agents, 
for instance neurons, that are characterized by a stochastic activity. It was hence 
important in that context to show that such orbits persist in these contexts, and 
to identify the mechanisms accounting for the emergence of such collective behav- 
iors. At the level of one single cell with stochastic perturbations, studies showed that 
systems presenting canard explosions displayed random MMOs, with the effect of ran- 
domly modifying the period of the oscillations and changing the qualitative nature 
of the solution even in small noise jB]. In this paper, the authors analyze stochas- 
tic perturbations of Fitzhugh-Nagumo model in the relaxation oscillations regime 
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(d) k = -2 (e) k = -3 



Figure 3.4. Solution to the mean-field equations \1.7\ (red and blue lines) and empirical 
averages of the network equation for N = 2000 neurons in the cases (a), (c), (d), (e) and (f) of 
Fig. \2.S\ In each case, we observe that qualitative behaviors are recovered in the network equations, 
and that early jumps occur systematically. 




concluded indeed noise resulted in sensible modifications of the trajectories with ran- 
dom interspike intervals and trajectories similar to MMOs with a random number 
of small oscillations that asymptotically distribute as a geometric random variable. 
As of mixed-mode oscillations, Berglund and collaborators [5j analyzed the effect of 
Gaussian noise on a 3 dimensional system with different timescales known to feature 
mixed-mode oscillations and show that noise changes the global dynamics and mixed 
mode patterns, and that there exists a maximal amount of noise above which small 
oscillations are hidden by noise. These studies pointed towards the fact that other 
phenomena, more resilient to noise, took place in large noisy systems giving rise to 



20 



M. DESROCHES, M. KRUPA & J. TOUBOUL 



such complex oscillations in the presence of multiple timescales. 

Here, we revisited these questions and proposed that self-averaging of the trajec- 
tories in large networks could provide a simple generic explanation to the emergence of 
complex oscillations in large-scale noisy networks. We based our analysis on a classical 
model of neuronal network, the Wilson-Cowan system, which describes heuristically 
the dynamics of firing-rates of neurons. This system presented two advantages in 
our view. First an unsolved mathematical problem was to rigorously derive asymp- 
totic limits of these networks with nonlinear mean-field interactions as the number 
of neurons tends to infinity. The nature of is indeed very different to usual network 
equations analyzed in the domain of mean-field limits and propagation of chaos as de- 
veloped in gas dynamics and extended to neurosciences [26l [28] , since the interacting 
term is a nonlinear function of a mean-field term, and moreover that the neurons are 
not governed by a stochastic differential equation but rather by a differential equation 
with stochastic coefficients Q. This lead us to extend the theory of mean- field limit^ 
towards this new kind of equations, which we developed in a relatively general fashion 
in section [ij Applied to Wilson-Cowan equations, this result allowed us to exactly 
reduce the network equations to a multiple timescales ODE acting on the mean of the 
solution. Interestingly, these equations reduce to a modified Wilson-Cowan system 
where the noise level is embedded in the nonlinearities. 

The analysis of the mean equation developed in section[2]had at least two interests. 
First, the fact that the asymptotic equation on the mean of the solution is itself 
solution of a Wilson-Cowan type of equation allowed to go deeper in the understanding 
of these models. And indeed, the behavior of Wilson-Cowan systems in the presence 
of multiple timescales was still an open problem important to solve since the original 
Wilson-Cowan model introduced in 1972 in .31j had an intrinsic slow-fast structure 
that was never explored as such in the literature. Second, this analysis allowed to 
understand the qualitative role of noise in the dynamics. As was done in [291 127] . 
we demonstrated here again that the noise did not have a pure disturbing effect on 
the trajectories but rather a dramatic qualitative effect on the form of the solutions, 
governing transitions from stationary to periodic behaviors and complex oscillatory 
behaviors. 

Of course, the same study would hold in the case of the simpler Amari-Hopfield 
model analyzed in , where the interaction between neurons is the sum of the firing- 
rate seen as a sigmoidal transform of the voltage. These models do not take into 
account the nonlinearities in the synaptic integration taken into account in Wilson- 
Cowan systems. For such networks, the mean-field equations has Gaussian asymptotic 
solutions, whose mean and standard deviation satisfy a closed set of ordinary differ- 
ential equations [5^ . Our analysis can easily be extended to such systems, where the 
limit equations obtained can be shown to have Gaussian solution whose mean satisfy 
the set of ODEs: 



T, 



P 



3=1 



where Gap is similar to our effective sigmoid function 1.12| and embeds dependence 



on the noise level ctq. In the presence of multiple timescales, we expect to observe 
the same kind of phenomena, and a similar dependence on noise levels, and these 



■^The proof uses similar tools as usual analysis. 
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phenomena may be rigorously demonstrated in the same flavor as we did in the present 
paper. 

Such rigorous results as proved in the present manuscript provide a proof of con- 
cept of the role of noise in large-scale multiple timescales systems. The question that 
may arise is whether similar phenomena arise in networks of excitable elements. In 
that case, the rigorous reduction to ODEs is no more possible, and proving rigorous 
results amounts studying infinite-dimensional systems given by non-local partial dif- 
ferential equations (see e.g. [28J for a formulation of the problem in excitable systems). 
However, preliminary numerical explorations of the problem tend to show that the 
same type of phenomena arise in these systems. The deep analysis of these systems 
and the existence of canards and complex oscillations in relationship with noise levels 
is a direct perspective of this work we are currently analyzing. 

Another deep problem raised by the present analysis is the systematic character- 
ization of the finite-size effects in bistable regimes. Indeed, we have observed that 
when there is a unique attractor of the mean- field equations, the finite-sized networks 
approximated fairly the predicted behavior, whereas in the bistable regime, switches 
occur randomly between the two attractors, but with a surprising regularity. Using 
the theory of the fluctuations around the mean-field limits [251 115) , we will explore in 
depth the origin of this regularity in a forthcoming study. 

The results obtained in this paper provide further insight on the role of noise in 
the dynamics of neuronal networks. More than a pure perturbation of the trajectories, 
noise shows to govern the dynamics of the networks, and moreover, the presence of 
canard explosions further assess that this dependence is extremely sensitive. This 
shows that stochastic networks displaying complex oscillations may be related to very 
finely tuned noise levels. This fine characterization of the noise may lead to unravel 
precisely noise levels in large-scale neuronal networks. 
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